#include "mydft.h"



int DFT_naive(Mat_<double>& gRe, Mat_<double>& gIm, Mat_<double>& GRe, Mat_<double>& GIm){
	const double pi = 3.1415926535897932384626433832795;
	const int N = gRe.rows;
	const int M = gRe.cols;
	const double MN = M*N;

	double a;
	double ca;
	double sa;

	//Scorro l'immagine trasformata
	for(int u = 0; u < N ; ++u){
		for(int v = 0; v < M ; ++v){
		
			GRe(u,v) = 0;
			GIm(u,v) = 0;
			//Per ogni pixel dell'immagine trasformata
			//Scorro l'immagine da trasformare
			for(int r = 0; r < N ; ++r){
				for(int c = 0; c < M ; ++c){
					
					a = -2 * pi * ( (u*r/double(N))   +   (v*c/(double(M)))) ;
					ca = cos(a);
					sa = sin(a);
					GRe(u,v) += gRe.at<double>(r,c) * ca - gIm.at<double>(r,c) * sa;
					GIm(u,v) += gRe.at<double>(r,c) * sa + gIm.at<double>(r,c) * ca;
				}
			
			} 	
		}
	}
	return 0;
}
int IDFT_naive(Mat_<double>& gRe, Mat_<double>& gIm, Mat_<double>& GRe, Mat_<double>& GIm){
	const double pi = 3.1415926535897932384626433832795;
	const int N = gRe.rows;
	const int M = gRe.cols;
	const double MN = M*N;

	double a;
	double ca;
	double sa;

	//Scorro l'immagine trasformata
	for(int u = 0; u < N ; ++u){
		for(int v = 0; v < M ; ++v){
		
			GRe(u,v) = 0;
			GIm(u,v) = 0;
			//Per ogni pixel dell'immagine trasformata
			//Scorro l'immagine da trasformare
			for(int r = 0; r < N ; ++r){
				for(int c = 0; c < M ; ++c){
					
					a = 2 * pi * ( (u*r/double(N))   +   (v*c/(double(M)))) ;
					ca = cos(a);
					sa = sin(a);
					GRe(u,v) += gRe.at<double>(r,c) * ca - gIm.at<double>(r,c) * sa;
					GIm(u,v) += gRe.at<double>(r,c) * sa + gIm.at<double>(r,c) * ca;
				}
			
			} 	
		}
	}

	
      for(int r = 0 ; r < N ; ++r){
          for(int c = 0 ; c < M ; ++c){
              GRe(r,c)/= MN;
              GIm(r,c)/= MN;
            }
        }
    
	return 0;
}
int DFT_separated(Mat_<double>& gRe, Mat_<double>& gIm, Mat_<double>& GRe, Mat_<double>& GIm){
  const double pi = 3.1415926535897932384626433832795;
  const int N = gRe.rows;
  const int M = gRe.cols;
  const double MN = M*N;

  Mat_<double> tmpRe(N,M);
  Mat_<double> tmpIm(N,M);

  //Eseguo la trasformata
  double a;
  double ca;
  double sa;

  
  for( int u = 0 ; u < M ; ++u){
	for(int v = 0 ; v < N ; ++v){
	//Per ogni riga u dell'immagine trasformata
		
		tmpRe(u,v) = 0;
		tmpIm(u,v) = 0;
		
		for(int r = 0; r < N; ++r){
			a = -2 * pi * u * r / double(N) ;
			ca = cos(a);
			sa = sin(a);
			tmpRe(u,v) += gRe.at<double>(r,v) * ca - gIm.at<double>(r,v) * sa;
			tmpIm(u,v) += gRe.at<double>(r,v) * sa + gIm.at<double>(r,v) * ca;	  
	  }
	}
  }


  for(int u = 0 ; u < M ; ++u){
	  for(int v = 0 ; v < N ; ++v){
		
		GRe(u,v) = 0;
		GIm(u,v) = 0;
		
		for(int c = 0 ; c < M ; ++c){
			a = -2 * pi * v * c / double(M) ;
			ca = cos(a);
			sa = sin(a);
			GRe(u,v) += tmpRe(u,c) * ca - tmpIm(u,c) * sa;
			GIm(u,v) += tmpRe(u,c) * sa + tmpIm(u,c) * ca;	  


		}

	  }
  }

 return 0;

}
int IDFT_separated(Mat_<double>& gRe, Mat_<double>& gIm, Mat_<double>& GRe, Mat_<double>& GIm){
  const double pi = 3.1415926535897932384626433832795;
  const int N = gRe.rows;
  const int M = gRe.cols;
  const double MN = M*N;

  Mat_<double> tmpRe(N,M);
  Mat_<double> tmpIm(N,M);

  //Eseguo la trasformata
  double a;
  double ca;
  double sa;

  
  for( int u = 0 ; u < M ; ++u){
	for(int v = 0 ; v < N ; ++v){
	//Per ogni riga u dell'immagine trasformata
		
		tmpRe(u,v) = 0;
		tmpIm(u,v) = 0;
		
		for(int r = 0; r < N; ++r){
			a = 2 * pi * u * r / double(N) ;
			ca = cos(a);
			sa = sin(a);
			tmpRe(u,v) += gRe.at<double>(r,v) * ca - gIm.at<double>(r,v) * sa;
			tmpIm(u,v) += gRe.at<double>(r,v) * sa + gIm.at<double>(r,v) * ca;	  
	  }
	}
  }


  for(int u = 0 ; u < M ; ++u){
	  for(int v = 0 ; v < N ; ++v){
		
		GRe(u,v) = 0;
		GIm(u,v) = 0;
		
		for(int c = 0 ; c < M ; ++c){
			a = 2 * pi * v * c / double(M) ;
			ca = cos(a);
			sa = sin(a);
			GRe(u,v) += tmpRe(u,c) * ca - tmpIm(u,c) * sa;
			GIm(u,v) += tmpRe(u,c) * sa + tmpIm(u,c) * ca;	  


		}

	  }
  }

 
      for(int u = 0 ; u < N ; ++u){
          for(int v = 0 ; v < M ; ++v){
              GRe(u,v)/= MN;
              GIm(u,v)/= MN;
            }
        }
    

  return 0;

}
int DFT_ptr(double* gRe ,double* gIm,double* GRe,double* GIm , int rows , int cols){
	
	//Costanti
	const double pi = 3.1415926535897932384626433832795;
	const int N = rows;
	const int M = cols;
	const int MN = rows*cols;

	//Immgini parziali
	double *tmpRe, *tmpIm;
	tmpRe = new double[MN];
	tmpIm = new double[MN];

	//Eseguo la trasformata
	register double a;
	double ca;
	double sa;



  //Calcolo per colonne
  for( int u = 0 ; u < M ; ++u){
	for(int v = 0 ; v < N ; ++v){
          tmpRe[u*M+v] = tmpIm[u*M+v] = 0;
          for(int r = 0; r < N ; ++r){
              a = - (2 * pi * u * r / double(N));
              ca = cos(a);
              sa = sin(a);
              tmpRe[u*M+v] += gRe[r*M+v] * ca - gIm[r*M+v] * sa;
              tmpIm[u*M+v] += gRe[r*M+v] * sa + gIm[r*M+v] * ca;
            }
        }
  }



  //Calcolo per righe
  for(int u = 0 ; u < M ; ++u){
	  for(int v = 0 ; v < N ; ++v){
          GRe[u*M+v] = gIm[u*M+v] = 0;
          for(int c = 0 ; c < M ; ++c){
             a = - (2 * pi * v * c / double(M));
             ca = cos(a);
             sa = sin(a);
             GRe[u*M+v] += tmpRe[u*M+c] * ca - tmpIm[u*M+c] * sa;
             GIm[u*M+v] += tmpRe[u*M+c] * sa + tmpIm[u*M+c] * ca;
            }
        }
    }
   return 0;
}
int IDFT_ptr(double* gRe ,double* gIm,double* GRe,double* GIm , int rows , int cols){
	
	//Costanti
	const double pi = 3.1415926535897932384626433832795;
	const int N = rows;
	const int M = cols;
	const int MN = rows*cols;

	//Immgini parziali
	double *tmpRe, *tmpIm;
	tmpRe = new double[MN];
	tmpIm = new double[MN];

	//Eseguo la trasformata
	register double a;
	double ca;
	double sa;



  //Calcolo per colonne
  for( int u = 0 ; u < M ; ++u){
	for(int v = 0 ; v < N ; ++v){
          tmpRe[u*M+v] = tmpIm[u*M+v] = 0;
          for(int r = 0; r < N ; ++r){
              a = 2 * pi * u * r / double(N);
              ca = cos(a);
              sa = sin(a);
              tmpRe[u*M+v] += gRe[r*M+v] * ca - gIm[r*M+v] * sa;
              tmpIm[u*M+v] += gRe[r*M+v] * sa + gIm[r*M+v] * ca;
            }
        }
  }

  //Calcolo per righe
  for(int u = 0 ; u < M ; ++u){
	  for(int v = 0 ; v < N ; ++v){
          GRe[u*M+v] = gIm[u*M+v] = 0;
          for(int c = 0 ; c < M ; ++c){
             a = 2 * pi * v * c / double(M);
             ca = cos(a);
             sa = sin(a);
             GRe[u*M+v] += tmpRe[u*M+c] * ca - tmpIm[u*M+c] * sa;
             GIm[u*M+v] += tmpRe[u*M+c] * sa + tmpIm[u*M+c] * ca;		 
            }
        }
    }
 
  for(int u = 0 ; u < M ; ++u){
	  for(int v = 0 ; v < N ; ++v){
		GRe[u*M+v]/= MN;
		GIm[u*M+v]/= MN;
	  }
  }
   
  return 0;
}

int FFT_separated(Mat_<float>& gRe, Mat_<float>& gIm, Mat_<float>& GRe, Mat_<float>& GIm){

  //Inizializzazione
	int N = gRe.rows;
	int M = gRe.cols;
	//Calcolo il log_2(N) e log_2(M) 
	int l2M = 0, p = 1; //l2M will become log_2(M)
	while(p < M) {p *= 2; l2M++;}
	int l2N = 0; p = 1; //l2N will become log_2(N)
	while(p < N) {p *= 2; l2N++;}
	
	M = 1 << l2M; N = 1 << l2N; //M,N devono essere potenze di 2 per funzionare
   
  //Inizializzo le immagini trasformate
  for(int r = 0 ; r < N ; ++r) {
	for(int c = 0 ; c < M ; ++c){
		GRe(r,c) = gRe(r,c);
		GIm(r,c) = gIm(r,c);
	}
  }

  //gRe.copyTo(GRe);
  //gIm.copyTo(GIm);


  //Eseguo l'inversione delle righe
  int j,k;
  for(int r = 0; r < N; ++r){ //per ogni riga
    j = 0;
    for(int i = 0; i < M - 1; ++i){
	  GRe(r,i) = gRe(r,j);
      GIm(r,i) = gIm(r,j);
      k = M / 2;
      while (k <= j) {j -= k; k/= 2;}
      j += k;
    }
  } 
  //Eseguo l'inversione delle colonne
  float tx = 0, ty = 0;  
  for(int c = 0; c < M; c++){ //per ogni colonna
	  j = 0;
    for(int i = 0; i < N - 1; i++){
      if(i < j){
        tx = GRe(i,c);
        ty = GIm(i,c);
        GRe(i,c) = GRe(j,c);
        GIm(i,c) = GIm(j,c);
        GRe(j,c) = tx;
        GIm(j,c) = ty;
      }
      k = N / 2;
      while (k <= j) {j -= k; k/= 2;}
      j += k;
    }
  }
  

  //Calcolo la FFT delle colonne
  float ca,sa,u1,u2,t1,t2,z;
  int l1,l2,i1;
  
  for(int c = 0; c < M; c++){ //per ogni colonna
    //FFT ad una dimensione:
    ca = -1.0;
    sa = 0.0;
    l1 = 1;
	l2 = 1;
    for(int l=0;l < l2M;l++){
      l1 = l2;
      l2 *= 2;
      u1 = 1.0;
      u2 = 0.0;
      for(int j = 0; j < l1; j++){
        for(int i = j; i < M; i += l2){
          i1 = i + l1;
          t1 = u1 * GRe(i1,c) - u2 * GIm(i1,c);
          t2 = u1 * GIm(i1,c) + u2 * GRe(i1,c);
          GRe(i1,c) = GRe(i,c) - t1;
          GIm(i1,c) = GIm(i,c) - t2;
          GRe(i,c) += t1;
          GIm(i,c) += t2;
        }
        z =  u1 * ca - u2 * sa;
        u2 = u1 * sa + u2 * ca;
        u1 = z;
      }
      sa = -sqrt((1.0 - ca) / 2.0);
      ca = sqrt((1.0 + ca) / 2.0);
    }
  }

 //Calcolo la FFT delle righe
  for(int r = 0; r < N; r++){ //per ogni riga
    //FFT ad una dimensione come prima:
    ca = -1.0;
    sa = 0.0;
    l1= 1;
	l2 = 1;
    for(int l = 0; l < l2N; l++){
      l1 = l2;
      l2 *= 2;
      u1 = 1.0;
      u2 = 0.0;
      for(int j = 0; j < l1; j++){
        for(int i = j; i < M; i += l2){
          int i1 = i + l1;
          t1 = u1 * GRe(r,i1) - u2 * GIm(r,i1);
          t2 = u1 * GIm(r,i1) + u2 * GRe(r,i1);
          GRe(r,i1) = GRe(r,i) - t1;
          GIm(r,i1) = GIm(r,i) - t2;
          GRe(r,i) += t1;
          GIm(r,i) += t2;
        }
        z =  u1 * ca - u2 * sa;
        u2 = u1 * sa + u2 * ca;
        u1 = z;
      }
      sa = -sqrt((1.0 - ca) / 2.0);
      ca = sqrt((1.0 + ca) / 2.0);
    }
  }
    
  //the values are divided through n if it isn't the inverse DFT.
  /*int d = N;
  for(int r = 0; r < N; r++) {
	  for(int c = 0; c < M; c++){
		GRe(r,c) /= d;
		GIm(r,c) /= d;
	  }
  }

  */

  return 0;

}
int IFFT_separated(Mat_<float>& gRe, Mat_<float>& gIm, Mat_<float>& GRe, Mat_<float>& GIm){

  //Inizializzazione
	int N = gRe.rows;
	int M = gRe.cols;
	//Calcolo il log_2(N) e log_2(M) 
	int l2M = 0, p = 1; //l2M will become log_2(M)
	while(p < M) {p *= 2; l2M++;}
	int l2N = 0; p = 1; //l2N will become log_2(N)
	while(p < N) {p *= 2; l2N++;}
	
	M = 1 << l2M; N = 1 << l2N; //M,N devono essere potenze di 2 per funzionare
   
  //Inizializzo le immagini trasformate
 // for(int r = 0 ; r < N ; ++r) {
	//for(int c = 0 ; c < M ; ++c){
	//	GRe(r,c) = gRe(r,c);
	//	GIm(r,c) = gIm(r,c);
	//}
 // }
  gRe.copyTo(GRe);
  gIm.copyTo(GIm);

  //Eseguo l'inversione delle righe
  int j,k;
  for(int r = 0; r < N; ++r){ //per ogni riga
    j = 0;
    for(int i = 0; i < M - 1; ++i){
      GRe(r,i) = gRe(r,j);
      GIm(r,i) = gIm(r,j);
      k = M / 2;
      while (k <= j) {j -= k; k/= 2;}
      j += k;
    }
  } 
  //Eseguo l'inversione delle colonne
  float tx = 0, ty = 0;  
  for(int c = 0; c < M; c++){ //per ogni colonna
	  j = 0;
    for(int i = 0; i < N - 1; i++){
      if(i < j){
        tx = GRe(i,c);
        ty = GIm(i,c);
        GRe(i,c) = GRe(j,c);
        GIm(i,c) = GIm(j,c);
        GRe(j,c) = tx;
        GIm(j,c) = ty;
      }
      k = N / 2;
      while (k <= j) {j -= k; k/= 2;}
      j += k;
    }
  }
  

  //Calcolo la FFT delle colonne
  float ca,sa,u1,u2,t1,t2,z;
  int l1,l2,i1;
  
  for(int c = 0; c < M; c++){ //per ogni colonna
    //FFT ad una dimensione:
    ca = -1.0;
    sa = 0.0;
    l1 = 1;
	l2 = 1;
    for(int l=0;l < l2M;l++){
      l1 = l2;
      l2 *= 2;
      u1 = 1.0;
      u2 = 0.0;
      for(int j = 0; j < l1; j++){
        for(int i = j; i < M; i += l2){
          i1 = i + l1;
          t1 = u1 * GRe(i1,c) - u2 * GIm(i1,c);
          t2 = u1 * GIm(i1,c) + u2 * GRe(i1,c);
          GRe(i1,c) = GRe(i,c) - t1;
          GIm(i1,c) = GIm(i,c) - t2;
          GRe(i,c) += t1;
          GIm(i,c) += t2;
        }
        z =  u1 * ca - u2 * sa;
        u2 = u1 * sa + u2 * ca;
        u1 = z;
      }
      sa = sqrt((1.0 - ca) / 2.0);
      ca = sqrt((1.0 + ca) / 2.0);
    }
  }

 //Calcolo la FFT delle righe
  for(int r = 0; r < N; r++){ //per ogni riga
    //FFT ad una dimensione come prima:
    ca = -1.0;
    sa = 0.0;
    l1= 1;
	l2 = 1;
    for(int l = 0; l < l2N; l++){
      l1 = l2;
      l2 *= 2;
      u1 = 1.0;
      u2 = 0.0;
      for(int j = 0; j < l1; j++){
        for(int i = j; i < M; i += l2){
          int i1 = i + l1;
          t1 = u1 * GRe(r,i1) - u2 * GIm(r,i1);
          t2 = u1 * GIm(r,i1) + u2 * GRe(r,i1);
          GRe(r,i1) = GRe(r,i) - t1;
          GIm(r,i1) = GIm(r,i) - t2;
          GRe(r,i) += t1;
          GIm(r,i) += t2;
        }
        z =  u1 * ca - u2 * sa;
        u2 = u1 * sa + u2 * ca;
        u1 = z;
      }
      sa = sqrt((1.0 - ca) / 2.0);
      ca = sqrt((1.0 + ca) / 2.0);
    }
  }
    
  //Divido ogni elemento per M*N
  int d = M*N;
  for(int r = 0; r < N; r++) {
	  for(int c = 0; c < M; c++){
		GRe(r,c) /= d;
		GIm(r,c) /= d;
	  }
  }
  
  return 0;

}
int FFT_ptr(float* gRe ,float* gIm, float* GRe, float* GIm , int rows , int cols){
	//Costanti
	int N = rows;
	int M = cols;

	//Calcolo il log_2(N) e log_2(M) 
	int l2M = 0, p = 1; 
	while(p < M) {p *= 2; l2M++;}
	int l2N = 0; p = 1;
	while(p < N) {p *= 2; l2N++;}
	
	M = 1 << l2M; N = 1 << l2N; //M,N devono essere potenze di 2 per funzionare
   
  //Inizializzo le immagini trasformate
  for(int r = 0 ; r < N ; ++r) {
	for(int c = 0 ; c < M ; ++c){
		GRe[r*M+c] = gRe[r*M+c];
		GIm[r*M+c] = gIm[r*M+c];
	}
  }
  
  //Eseguo l'inversione delle righe
  int j,k;
  for(int r = 0; r < N; ++r){ //per ogni riga
    j = 0;
    for(int i = 0; i < M - 1; ++i){
      GRe[r*M+i] = gRe[r*M+j];
      GIm[r*M+i] = gIm[r*M+j];
      k = M / 2;
      while (k <= j) {j -= k; k/= 2;}
      j += k;
    }
  } 
  //Eseguo l'inversione delle colonne
  float tx = 0, ty = 0;  
  for(int c = 0; c < M; c++){ //per ogni colonna
	  j = 0;
    for(int i = 0; i < N - 1; i++){
      if(i < j){
        tx = GRe[i*M+c];
        ty = GIm[i*M+c];
        GRe[i*M+c] = GRe[j*M+c];
        GIm[i*M+c] = GIm[j*M+c];
        GRe[j*M+c] = tx;
        GIm[j*M+c] = ty;
      }
      k = N / 2;
      while (k <= j) {j -= k; k/= 2;}
      j += k;
    }
  }
  
  

  //Calcolo la FFT delle colonne
  float ca,sa,u1,u2,t1,t2,z;
  int l1,l2,i1;
  
  for(int c = 0; c < M; c++){ //per ogni colonna
    //FFT ad una dimensione:
    ca = -1.0;
    sa = 0.0;
    l1 = 1;
	l2 = 1;
    for(int l=0;l < l2M;l++){
      l1 = l2;
      l2 *= 2;
      u1 = 1.0;
      u2 = 0.0;
      for(int j = 0; j < l1; j++){
        for(int i = j; i < M; i += l2){
          i1 = i + l1;
          t1 = u1 * GRe[i1*M+c] - u2 * GIm[i1*M+c];
          t2 = u1 * GIm[i1*M+c] + u2 * GRe[i1*M+c];
          GRe[i1*M+c] = GRe[i*M+c] - t1;
          GIm[i1*M+c] = GIm[i*M+c] - t2;
          GRe[i*M+c] += t1;
          GIm[i*M+c] += t2;
        }
        z =  u1 * ca - u2 * sa;
        u2 = u1 * sa + u2 * ca;
        u1 = z;
      }
      sa = -sqrt((1.0 - ca) / 2.0);
      ca = sqrt((1.0 + ca) / 2.0);
    }
  }

 //Calcolo la FFT delle righe
  for(int r = 0; r < N; r++){ //per ogni riga
    //FFT ad una dimensione come prima:
    ca = -1.0;
    sa = 0.0;
    l1= 1;
	l2 = 1;
    for(int l = 0; l < l2N; l++){
      l1 = l2;
      l2 *= 2;
      u1 = 1.0;
      u2 = 0.0;
      for(int j = 0; j < l1; j++){
        for(int i = j; i < M; i += l2){
          int i1 = i + l1;
          t1 = u1 * GRe[r*M+i1] - u2 * GIm[r*M+i1];
          t2 = u1 * GIm[r*M+i1] + u2 * GRe[r*M+i1];
          GRe[r*M+i1] = GRe[r*M+i] - t1;
          GIm[r*M+i1] = GIm[r*M+i] - t2;
          GRe[r*M+i] += t1;
          GIm[r*M+i] += t2;
        }
        z =  u1 * ca - u2 * sa;
        u2 = u1 * sa + u2 * ca;
        u1 = z;
      }
      sa = -sqrt((1.0 - ca) / 2.0);
      ca = sqrt((1.0 + ca) / 2.0);
    }
  }
    
  return 0;

}
int IFFT_ptr(float* gRe ,float* gIm, float* GRe, float* GIm , int rows , int cols){
	//Costanti
	int N = rows;
	int M = cols;

	//Calcolo il log_2(N) e log_2(M) 
	int l2M = 0, p = 1; 
	while(p < M) {p *= 2; l2M++;}
	int l2N = 0; p = 1;
	while(p < N) {p *= 2; l2N++;}
	
	M = 1 << l2M; N = 1 << l2N; //M,N devono essere potenze di 2 per funzionare
   
  //Inizializzo le immagini trasformate
  for(int r = 0 ; r < N ; ++r) {
	for(int c = 0 ; c < M ; ++c){
		GRe[r*M+c] = gRe[r*M+c];
		GIm[r*M+c] = gIm[r*M+c];
	}
  }

  //Eseguo l'inversione delle righe
  int j,k;
  for(int r = 0; r < N; ++r){ //per ogni riga
    j = 0;
    for(int i = 0; i < M - 1; ++i){
      GRe[r*M+i] = gRe[r*M+j];
      GIm[r*M+i] = gIm[r*M+j];
      k = M / 2;
      while (k <= j) {j -= k; k/= 2;}
      j += k;
    }
  } 
  //Eseguo l'inversione delle colonne
  float tx = 0, ty = 0;  
  for(int c = 0; c < M; c++){ //per ogni colonna
	  j = 0;
    for(int i = 0; i < N - 1; i++){
      if(i < j){
        tx = GRe[i*M+c];
        ty = GIm[i*M+c];
        GRe[i*M+c] = GRe[j*M+c];
        GIm[i*M+c] = GIm[j*M+c];
        GRe[j*M+c] = tx;
        GIm[j*M+c] = ty;
      }
      k = N / 2;
      while (k <= j) {j -= k; k/= 2;}
      j += k;
    }
  }
  

  //Calcolo la FFT delle colonne
  float ca,sa,u1,u2,t1,t2,z;
  int l1,l2,i1;
  
  for(int c = 0; c < M; c++){ //per ogni colonna
    //FFT ad una dimensione:
    ca = -1.0;
    sa = 0.0;
    l1 = 1;
	l2 = 1;
    for(int l=0;l < l2M;l++){
      l1 = l2;
      l2 *= 2;
      u1 = 1.0;
      u2 = 0.0;
      for(int j = 0; j < l1; j++){
        for(int i = j; i < M; i += l2){
          i1 = i + l1;
          t1 = u1 * GRe[i1*M+c] - u2 * GIm[i1*M+c];
          t2 = u1 * GIm[i1*M+c] + u2 * GRe[i1*M+c];
          GRe[i1*M+c] = GRe[i*M+c] - t1;
          GIm[i1*M+c] = GIm[i*M+c] - t2;
          GRe[i*M+c] += t1;
          GIm[i*M+c] += t2;
        }
        z =  u1 * ca - u2 * sa;
        u2 = u1 * sa + u2 * ca;
        u1 = z;
      }
      sa = sqrt((1.0 - ca) / 2.0);
      ca = sqrt((1.0 + ca) / 2.0);
    }
  }

 //Calcolo la FFT delle righe
  for(int r = 0; r < N; r++){ //per ogni riga
    //FFT ad una dimensione come prima:
    ca = -1.0;
    sa = 0.0;
    l1= 1;
	l2 = 1;
    for(int l = 0; l < l2N; l++){
      l1 = l2;
      l2 *= 2;
      u1 = 1.0;
      u2 = 0.0;
      for(int j = 0; j < l1; j++){
        for(int i = j; i < M; i += l2){
          int i1 = i + l1;
          t1 = u1 * GRe[r*M+i1] - u2 * GIm[r*M+i1];
          t2 = u1 * GIm[r*M+i1] + u2 * GRe[r*M+i1];
          GRe[r*M+i1] = GRe[r*M+i] - t1;
          GIm[r*M+i1] = GIm[r*M+i] - t2;
          GRe[r*M+i] += t1;
          GIm[r*M+i] += t2;
        }
        z =  u1 * ca - u2 * sa;
        u2 = u1 * sa + u2 * ca;
        u1 = z;
      }
      sa = sqrt((1.0 - ca) / 2.0);
      ca = sqrt((1.0 + ca) / 2.0);
    }
  }

  //Divido ogni elemento per M*N
  int d = M*N;
  for(int r = 0; r < N; r++) {
	  for(int c = 0; c < M; c++){
		GRe[r*M+c] /= d;
		GIm[r*M+c] /= d;
	  }
  }
    
  return 0;

}





/*

int testDFT(){
	
	Mat3b gg = imread("limbo1024.jpeg");
	
	Mat_<double> g(gg.rows,gg.cols);
	cvtColor(gg,g,CV_BGR2GRAY);
	namedWindow("gRe");
	imshow("gRe",g);
	//Creo immagini da elaborare
  Mat_<double> gRe(g.rows,g.cols);
  g.convertTo(gRe,CV_64F);
  Mat_<double> gIm(g.rows,g.cols);
  Mat_<double> GRe(g.rows,g.cols);
  Mat_<double> GIm(g.rows,g.cols);
  
  Mat_<double> gReNEW(g.rows,g.cols);
  Mat_<double> gImNEW(g.rows,g.cols);
  clock_t startTime;
  
  
   ////////////DFT Naive
  gIm = Mat::zeros(g.rows,g.cols,sizeof(double));
  GRe = Mat::zeros(g.rows,g.cols,sizeof(double));
  GIm = Mat::zeros(g.rows,g.cols,sizeof(double));
  startTime = clock();//
  DFT_naive(gRe,gIm,GRe,GIm);
  IDFT_naive(GRe,GIm,gReNEW,gImNEW);
  cout << "DFT Naive: "<< double( clock() - startTime ) / (double)CLOCKS_PER_SEC<< " secondi." << endl;

  ////////////DFT Separated
  gIm = Mat::zeros(g.rows,g.cols,sizeof(double));
  GRe = Mat::zeros(g.rows,g.cols,sizeof(double));
  GIm = Mat::zeros(g.rows,g.cols,sizeof(double));
  startTime = clock();
  DFT_separated(gRe,gIm,GRe,GIm);
  IDFT_separated(GRe,GIm,gReNEW,gImNEW);
  cout << "DFT Separated: "<< double( clock() - startTime ) / (double)CLOCKS_PER_SEC<< " secondi." << endl;
  
  ///////////DFT Puntatori
  gIm = Mat::zeros(g.rows,g.cols,sizeof(double));
  GRe = Mat::zeros(g.rows,g.cols,sizeof(double));
  GIm = Mat::zeros(g.rows,g.cols,sizeof(double));
  startTime = clock();
  DFT_ptr( (double*)gRe.data , (double*)gIm.data, (double*)GRe.data ,(double*)GIm.data ,gRe.rows, gRe.cols);
  IDFT_ptr( (double*)GRe.data , (double*)GIm.data, (double*)gReNEW.data ,(double*)gImNEW.data ,GRe.rows, GRe.cols);
  cout << "DFT Puntatori: "<< double( clock() - startTime ) / (double)CLOCKS_PER_SEC<< " secondi." << endl;
 //  * /
  ///////////FFT Separated
  gIm = Mat::zeros(g.rows,g.cols,sizeof(double));
  GRe = Mat::zeros(g.rows,g.cols,sizeof(double));
  GIm = Mat::zeros(g.rows,g.cols,sizeof(double));
  startTime = clock();
  FFT_separated(gRe,gIm,GRe,GIm);
  IFFT_separated(GRe,GIm,gReNEW,gImNEW);
  cout << "FFT Separated: "<< double( clock() - startTime ) / (double)CLOCKS_PER_SEC<< " secondi." << endl;

  ///////////FFT Puntatori
  gIm = Mat::zeros(g.rows,g.cols,sizeof(double));
  GRe = Mat::zeros(g.rows,g.cols,sizeof(double));
  GIm = Mat::zeros(g.rows,g.cols,sizeof(double));
  startTime = clock();
  FFT_ptr( (double*)gRe.data , (double*)gIm.data, (double*)GRe.data ,(double*)GIm.data ,gRe.rows, gRe.cols);
  IFFT_ptr( (double*)GRe.data , (double*)GIm.data, (double*)gReNEW.data ,(double*)gImNEW.data ,GRe.rows, GRe.cols);
  cout << "FFT Puntatori: "<< double( clock() - startTime ) / (double)CLOCKS_PER_SEC<< " secondi." << endl;
  
  Mat_<uchar> outRe(g.rows,g.cols);
  gReNEW.convertTo(outRe,CV_8U);
  namedWindow("outRe");
  imshow("outRe",outRe);

#ifdef DEBUG_ON
  //MOSTRO IMMAGINI TRASFORMATE
  Mat1b GReLOG(GRe.rows,GRe.cols);
  Mat1b GImLOG(GRe.rows,GRe.cols);
  Mat1b::iterator j;
  Mat_<double>::iterator i;
  for(i = GRe.begin(), j = GReLOG.begin(); i != GRe.end(); ++i, ++j)
      (*j) = log(double(1+abs(*i)));
  for(i = GIm.begin(), j = GImLOG.begin(); i != GIm.end(); ++i, ++j)
      (*j) = log(double(1+abs(*i)));
  namedWindow("GReLOG");
  imshow("GReLOG",GReLOG);
  namedWindow("GImLOG");
  imshow("GImLOG",GImLOG);
#endif
  
  
  waitKey();
  return 0;
}

*/